subroutine determine_multisnap(isnap,is_multi)
  use hdf5_wrapper
  use logs
  use myfiles
  implicit none
  integer(hid_t)       :: file_id,error  
  integer*4,intent(in) :: isnap
  character(len=16)    :: csnap
  character(len=512)   :: infile_no_ext,fname
  logical,intent(out)  :: is_multi

  ! Turn off hdf5 error messages
  call h5eset_auto_f(0,error)  
  
  write(csnap,'(I3.3)')isnap
  infile_no_ext = trim(adjustl(indir))//'/SelectedParticles_'//trim(adjustl(csnap))//'.hdf5'

  ! Test for presence of file - if not there, add on .0 to the filename
  fname = trim(adjustl(infile_no_ext))//'.hdf5'
  
  call h5fopen_f(fname,H5F_ACC_RDONLY_F,file_id,error)
  if(error.lt.0) then
     ! Can't open, try .0.hdf5 instead and try again
     fname = trim(adjustl(infile_no_ext))//'.0.hdf5'
     call h5fopen_f(fname,H5F_ACC_RDONLY_F,file_id,error)
     if(error < 0) then
        call abort('cannot find snapshot in single- or multi-file format: '//trim(adjustl(infile_no_ext)))
     else  
        is_multi = .TRUE.
        call h5fclose_f(file_id,error)
     endif
  else     
     is_multi = .FALSE.
     call h5fclose_f(file_id,error)
  endif
  
  ! Turn on hdf5 error messages
  call h5eset_auto_f(1,error)
  
end subroutine determine_multisnap
          
  
subroutine read_snapshot_header(isnap,ifile)
  use hdf5_wrapper
  use header_vars
  use myfiles
  implicit none
  integer             :: snap_handle
  integer*4           :: isnap,ifile
  character(len=16)   :: csnap,cfile
  character(len=512)  :: infile_no_ext,fname
  logical             :: is_multi


  call determine_multisnap(isnap,is_multi)

  write(csnap,'(I3.3)')isnap
  infile_no_ext = trim(adjustl(indir))//'/SelectedParticles_'//trim(adjustl(csnap))//'.hdf5'

  if(is_multi) then 
     write(cfile,'(I)')ifile
     fname = trim(adjustl(infile_no_ext))//'.'//trim(adjustl(cfile))//'.hdf5'
  else
     fname = trim(adjustl(infile_no_ext))//'.hdf5'
  endif

  call hdf5_open_file(snap_handle,trim(fname),readonly=.true.)
  call hdf5_read_attribute(snap_handle,'Header/NumPart_ThisFile',header%NumPart_ThisFile)
  call hdf5_read_attribute(snap_handle,'Header/NumPart_Total',header%NumPart_Total)
  call hdf5_read_attribute(snap_handle,'Header/NumPart_Total_HighWord',header%NumPart_Total_HighWord)
  call hdf5_read_attribute(snap_handle,'Header/MassTable',header%MassTable)
  call hdf5_read_attribute(snap_handle,'Header/ExpansionFactor',header%ExpansionFactor)
  call hdf5_read_attribute(snap_handle,'Header/Time_GYR',header%Time_GYR)
  call hdf5_read_attribute(snap_handle,'Header/Redshift',header%Redshift)
  call hdf5_read_attribute(snap_handle,'Header/BoxSize',header%BoxSize)
  call hdf5_read_attribute(snap_handle,'Header/NumFilesPerSnapshot',header%NumFilesPerSnapshot)
  call hdf5_read_attribute(snap_handle,'Header/Omega0',header%Omega0)
  call hdf5_read_attribute(snap_handle,'Header/OmegaLambda',header%OmegaLambda)
  call hdf5_read_attribute(snap_handle,'Header/OmegaBaryon',header%OmegaBaryon)
  call hdf5_read_attribute(snap_handle,'Header/HubbleParam',header%HubbleParam)
  call hdf5_read_attribute(snap_handle,'Header/Flag_Sfr',header%Flag_Sfr)
  call hdf5_read_attribute(snap_handle,'Header/Flag_Cooling',header%Flag_Cooling)
  call hdf5_read_attribute(snap_handle,'Header/Flag_StellarAge',header%Flag_StellarAge)
  call hdf5_read_attribute(snap_handle,'Header/Flag_Metals',header%Flag_Metals)
  call hdf5_read_attribute(snap_handle,'Header/Flag_Feedback',header%Flag_Feedback)
  call hdf5_close_file(snap_handle)

end subroutine read_snapshot_header

subroutine read_original_selection_header()
  use hdf5_wrapper
  use myfiles
  use particle_data
  use param_module
  implicit none
  integer             :: snap_handle
  integer*4           :: isnap
  character(len=16)   :: csnap
  character(len=512)  :: fname
 
  write(csnap,'(I3.3)')iend
  fname = trim(adjustl(indir))//'/OriginalSelection_'//trim(adjustl(csnap))//'.hdf5'

  call hdf5_open_file(snap_handle,fname, readonly=.true.)
  call hdf5_read_attribute(snap_handle,'Header/NumSelected',NumSelected)
  call hdf5_read_attribute(snap_handle,'Header/x0'         ,x0)
  call hdf5_read_attribute(snap_handle,'Header/y0'         ,y0)
  call hdf5_read_attribute(snap_handle,'Header/z0'         ,z0)
  call hdf5_read_attribute(snap_handle,'Header/r0'         ,r0)
  call hdf5_close_file(snap_handle)

end subroutine read_original_selection_header

subroutine read_select_file_header(isnap)
  use hdf5_wrapper
  use header_vars
  use myfiles
  implicit none
  integer             :: snap_handle
  integer*4           :: isnap
  character(len=16)   :: csnap
  character(len=512)  :: fname
 
  write(csnap,'(I3.3)')isnap
  fname = trim(adjustl(indir))//'/SelectedParticles_'//trim(adjustl(csnap))//'.hdf5'

  call hdf5_open_file(snap_handle,fname, readonly=.true.)
  call hdf5_read_attribute(snap_handle,'Header/NumGathered',NumGathered)
  call hdf5_close_file(snap_handle)

  ! SPOOF
  !NumGathered = 0

end subroutine read_select_file_header

subroutine read_all_ids()
  use logs
  use hdf5_wrapper
  use header_vars
  use param_module
  use myfiles
  use particle_data

  implicit none
  integer                   :: snap_handle
  integer*4                 :: isnap,ifile,itype,i,j
  integer*4, dimension(0:5) :: Offset = 0
  character(len=16)         :: csnap,cfile
  character(len=512)        :: infile_no_ext,fname
  logical                   :: is_multi  

  write(logunit,*)' -reading IDs from OriginalSelection...'
  call flush(logunit)

  write(csnap,'(I3.3)')iend
  fname = trim(adjustl(indir))//'/OriginalSelection_'//trim(adjustl(csnap))//'.hdf5'

  Offset(:) = 0

  call hdf5_open_file(snap_handle,fname,readonly=.true.)
  i = Offset(0)
  j = Offset(0)+NumSelected(0)-1
  if(NumSelected(0) > 0) call hdf5_read_data(snap_handle,'PartType0/ParticleIDs',ParticleIDs%Type0(i:j))
  if(NumSelected(0) > 0) SnapIDs%Type0(i:j) = iend
  i = Offset(1)
  j = Offset(1)+NumSelected(1)-1
  if(NumSelected(1) > 0) call hdf5_read_data(snap_handle,'PartType1/ParticleIDs',ParticleIDs%Type1(i:j))
  if(NumSelected(1) > 0) SnapIDs%Type1(i:j) = iend
  i = Offset(2)
  j = Offset(2)+NumSelected(2)-1
  if(NumSelected(2) > 0) call hdf5_read_data(snap_handle,'PartType2/ParticleIDs',ParticleIDs%Type2(i:j))
  if(NumSelected(2) > 0) SnapIDs%Type2(i:j) = iend
  i = Offset(3)
  j = Offset(3)+NumSelected(3)-1
  if(NumSelected(3) > 0) call hdf5_read_data(snap_handle,'PartType3/ParticleIDs',ParticleIDs%Type3(i:j))
  if(NumSelected(3) > 0) SnapIDs%Type3(i:j) = iend
  i = Offset(4)
  j = Offset(4)+NumSelected(4)-1
  if(NumSelected(4) > 0) call hdf5_read_data(snap_handle,'PartType4/ParticleIDs',ParticleIDs%Type4(i:j))
  if(NumSelected(4) > 0) SnapIDs%Type4(i:j) = iend
  i = Offset(5)
  j = Offset(5)+NumSelected(5)-1
  if(NumSelected(5) > 0) call hdf5_read_data(snap_handle,'PartType5/ParticleIDs',ParticleIDs%Type5(i:j))
  if(NumSelected(5) > 0) SnapIDs%Type5(i:j) = iend
   call hdf5_close_file(snap_handle)
   
   Offset(0:5) = Offset(0:5) + NumSelected(0:5)

  write(logunit,*)' -reading IDs from SelectedParticles...'
  call flush(logunit)

  ! Loop over all snapshots, all files, and all types
  do isnap = istart,iend

     write(csnap,'(I3.3)')isnap
     fname = trim(adjustl(indir))//'/SelectedParticles_'//trim(adjustl(csnap))//'.hdf5'

     call read_select_file_header(isnap)
     call hdf5_open_file(snap_handle,fname,readonly=.true.)
     i = Offset(0)
     j = Offset(0)+NumGathered(0)-1
     if(NumGathered(0) > 0) call hdf5_read_data(snap_handle,'PartType0/ParticleIDs',ParticleIDs%Type0(i:j))
     if(NumGathered(0) > 0) SnapIDs%Type0(i:j) = isnap
     i = Offset(1)
     j = Offset(1)+NumGathered(1)-1
     if(NumGathered(1) > 0) call hdf5_read_data(snap_handle,'PartType1/ParticleIDs',ParticleIDs%Type1(i:j))
     if(NumGathered(1) > 0) SnapIDs%Type1(i:j) = isnap
     i = Offset(2)
     j = Offset(2)+NumGathered(2)-1
     if(NumGathered(2) > 0) call hdf5_read_data(snap_handle,'PartType2/ParticleIDs',ParticleIDs%Type2(i:j))
     if(NumGathered(2) > 0) SnapIDs%Type2(i:j) = isnap
     i = Offset(3)
     j = Offset(3)+NumGathered(3)-1
     if(NumGathered(3) > 0) call hdf5_read_data(snap_handle,'PartType3/ParticleIDs',ParticleIDs%Type3(i:j))
     if(NumGathered(3) > 0) SnapIDs%Type3(i:j) = isnap
     i = Offset(4)
     j = Offset(4)+NumGathered(4)-1
     if(NumGathered(4) > 0) call hdf5_read_data(snap_handle,'PartType4/ParticleIDs',ParticleIDs%Type4(i:j))
     if(NumGathered(4) > 0) SnapIDs%Type4(i:j) = isnap
     i = Offset(5)
     j = Offset(5)+NumGathered(5)-1
     if(NumGathered(5) > 0) call hdf5_read_data(snap_handle,'PartType5/ParticleIDs',ParticleIDs%Type5(i:j))
     if(NumGathered(5) > 0) SnapIDs%Type5(i:j) = isnap
     call hdf5_close_file(snap_handle)

     Offset(0:5) = Offset(0:5) + NumGathered(0:5)
     
  enddo

end subroutine read_all_ids
       

subroutine read_binary_ics()
  use logs
  use hdf5_wrapper
  use header_vars
  use param_module
  use memory_module
  use myfiles
  use particle_data
  use ic_data

  implicit none
  integer                              :: snap_handle,myunit
  integer*4                            :: isnap,ifile,itype
  integer*4                            :: i,j,itemp,jtemp,idim
  integer*4                            :: ios,ntemp
  integer*4, dimension(0:5)            :: Offset = 0
  integer*4                            :: TempOffset = 0
  character(len=16)                    :: csnap,cfile
  character(len=512)                   :: infile_no_ext,fname
  logical                              :: is_multi  
  logical                              :: byteswap
  real*4, allocatable, dimension(:,:)  :: pos
  real*4, allocatable, dimension(:,:)  :: vel
  integer*4, allocatable, dimension(:) :: ids

  
  write(logunit,*)' -reading binary initial conditions'
  call flush(logunit)

  myunit = 50

  ! Try opening single file
  infile_no_ext = trim(adjustl(icdir))//trim(adjustl(icfile))
  fname = trim(adjustl(infile_no_ext))
  open(myunit,file=fname,form='unformatted',status='old',iostat=ios,access='stream')
  if(ios /= 0)then 
     ! Try first of multifile
     fname = trim(adjustl(infile_no_ext))//'.0'
     open(myunit,file=fname,form='unformatted',status='old',iostat=ios,access='stream')
     if(ios /= 0) call abort(" -cannot read initial conditions in binary format")
     is_multi = .TRUE.
     read(myunit)ntemp
     close(myunit)
  else
     is_multi = .FALSE.
     read(myunit)ntemp
     close(myunit)
  endif

  ! Might we need to byteswap?
  write(logunit,*)'  -nbytes(header) reported as:',ntemp
  if(ntemp /= 256)then
     write(logunit,*)'  -will byteswap'
     byteswap = .TRUE.
     myunit   = 51
  else
     write(logunit,*)'  -no need to byteswap'
     byteswap = .FALSE.
  endif

  fname = trim(adjustl(infile_no_ext))
  if(is_multi) fname = trim(adjustl(infile_no_ext))//'.0'

  open(myunit,file=fname,form='unformatted',status='old',iostat=ios)
  read(myunit)binary_ics_header
  ! Allocate memory for particle data
  NumPart_ICs = binary_ics_header%NumPart_Total
  write(logunit,*)' -Initial Conditions comprise ',binary_ics_header%NumFilesPerSnapshot,' files and:'
  write(logunit,'(A11,I9.9)')'  -Type0 = ',NumPart_ICs(0)
  write(logunit,'(A11,I9.9)')'  -Type1 = ',NumPart_ICs(1)
  write(logunit,'(A11,I9.9)')'  -Type2 = ',NumPart_ICs(2)
  write(logunit,'(A11,I9.9)')'  -Type3 = ',NumPart_ICs(3)
  write(logunit,'(A11,I9.9)')'  -Type4 = ',NumPart_ICs(4)
  write(logunit,'(A11,I9.9)')'  -Type5 = ',NumPart_ICs(5)

  call allocate_ics_particles()
  close(myunit)

  ! Now loop over files
  do ifile=0,binary_ics_header%NumFilesPerSnapshot-1

     ! Report progress
     write(cfile,'(I)')ifile
     if(ifile /= binary_ics_header%NumFilesPerSnapshot-1)then
        write (logunit,'(A)',advance='no')char(13)//'  -file = '//trim(adjustl(cfile))
        call flush(logunit)
     else
        write (logunit,'(A)',advance='yes')char(13)//'  -file = '//trim(adjustl(cfile))
        call flush(logunit)
     endif

     fname = trim(adjustl(infile_no_ext))
     if(is_multi) fname = trim(adjustl(infile_no_ext))//'.'//trim(adjustl(cfile))

     open(myunit,file=fname,form='unformatted',status='old',iostat=ios)
     read(myunit)binary_ics_header

     if(allocated(pos)) deallocate(pos)
     allocate(pos(0:2,0:sum(binary_ics_header%NumPart_ThisFile)-1))
     if(.not.(allocated(pos))) call abort(" -failed to allocate pos")

     if(allocated(vel)) deallocate(vel)
     allocate(vel(0:2,0:sum(binary_ics_header%NumPart_ThisFile)-1))
     if(.not.(allocated(vel))) call abort(" -failed to allocate vel")

     if(allocated(ids)) deallocate(ids)
     allocate(ids(0:sum(binary_ics_header%NumPart_ThisFile)-1))
     if(.not.(allocated(ids))) call abort(" -failed to allocate ids")

     ! Coordinates
     read(myunit)pos
     read(myunit)vel
     read(myunit)ids

     itype = 0
     i     = Offset(itype)
     j     = Offset(itype)  + binary_ics_header%NumPart_ThisFile(itype)-1
     itemp = 0
     jtemp = itemp + binary_ics_header%NumPart_ThisFile(itype)-1
     do idim=0,2 
        IC_Coordinates%Type0(idim,i:j) = pos(idim,itemp:jtemp)
     enddo
     do idim=0,2 
        IC_Velocity%Type0(idim,i:j)    = vel(idim,itemp:jtemp)
     enddo
     IC_ParticleIDs%Type0(i:j)     = ids(itemp:jtemp)

     itype = 1
     i     = Offset(itype)
     j     = Offset(itype)  + binary_ics_header%NumPart_ThisFile(itype)-1
     itemp = jtemp + 1
     jtemp = itemp + binary_ics_header%NumPart_ThisFile(itype)-1
     do idim=0,2 
        IC_Coordinates%Type1(idim,i:j) = pos(idim,itemp:jtemp)
     enddo
     do idim=0,2 
        IC_Velocity%Type1(idim,i:j)    = vel(idim,itemp:jtemp)
     enddo
     IC_ParticleIDs%Type1(i:j)     = ids(itemp:jtemp)

     itype = 2
     i     = Offset(itype)
     j     = Offset(itype)  + binary_ics_header%NumPart_ThisFile(itype)-1
     itemp = jtemp + 1
     jtemp = itemp + binary_ics_header%NumPart_ThisFile(itype)-1
     do idim=0,2
        IC_Coordinates%Type2(idim,i:j) = pos(idim,itemp:jtemp)
     enddo
     do idim=0,2
        IC_Velocity%Type2(idim,i:j)    = vel(idim,itemp:jtemp)
     enddo
     IC_ParticleIDs%Type2(i:j)     = ids(itemp:jtemp)

     itype = 3
     i     = Offset(itype)
     j     = Offset(itype)  + binary_ics_header%NumPart_ThisFile(itype)-1
     itemp = jtemp + 1
     jtemp = itemp + binary_ics_header%NumPart_ThisFile(itype)-1
     do idim=0,2
        IC_Coordinates%Type3(idim,i:j) = pos(idim,itemp:jtemp)
     enddo
     do idim=0,2
        IC_Velocity%Type3(idim,i:j)    = vel(idim,itemp:jtemp)
     enddo
     IC_ParticleIDs%Type3(i:j)     = ids(itemp:jtemp)

     itype = 4
     i     = Offset(itype)
     j     = Offset(itype)  + binary_ics_header%NumPart_ThisFile(itype)-1
     itemp = jtemp + 1
     jtemp = itemp + binary_ics_header%NumPart_ThisFile(itype)-1
     do idim=0,2
        IC_Coordinates%Type4(idim,i:j) = pos(idim,itemp:jtemp)
     enddo
     do idim=0,2
        IC_Velocity%Type4(idim,i:j)    = vel(idim,itemp:jtemp)
     enddo
     IC_ParticleIDs%Type4(i:j)     = ids(itemp:jtemp)

     itype = 5
     i     = Offset(itype)
     j     = Offset(itype)  + binary_ics_header%NumPart_ThisFile(itype)-1
     itemp = jtemp + 1
     jtemp = itemp + binary_ics_header%NumPart_ThisFile(itype)-1
     do idim=0,2
        IC_Coordinates%Type5(idim,i:j) = pos(idim,itemp:jtemp)
     enddo
     do idim=0,2
        IC_Velocity%Type5(idim,i:j)    = vel(idim,itemp:jtemp)
     enddo
     IC_ParticleIDs%Type5(i:j)     = ids(itemp:jtemp)

     ! Update global offsets
     do itype = 0,5
        Offset(itype) = Offset(itype) + binary_ics_header%NumPart_ThisFile(itype)
     enddo

  enddo

end subroutine read_binary_ics




subroutine read_hdf5_ics()
  use logs
  use hdf5_wrapper
  use header_vars
  use param_module
  use myfiles
  use particle_data

  implicit none
  integer                   :: snap_handle
  integer*4                 :: isnap,ifile,itype,i,j
  integer*4, dimension(0:5) :: Offset,FileOffset = 0
  character(len=16)         :: csnap,cfile
  character(len=512)        :: infile_no_ext,fname
  logical                   :: is_multi  


  write(logunit,*)' -reading hdf5 initial conditions'
  call flush(logunit)

  ! Loop over all snapshots, all files, and all types
  do isnap = istart,iend

     ifile = 0
     call read_snapshot_header(isnap,ifile)
     call determine_multisnap(isnap,is_multi)

     write(csnap,'(I3.3)')isnap
     infile_no_ext = trim(adjustl(indir))//'/Selected.snap_'//trim(adjustl(csnap))

     do itype = 0,5
        Offset(itype) = SnapOffset(itype,isnap-istart)
     enddo

     do ifile = 0,header%NumFilesPerSnapshot-1

        if(is_multi) then 
           write(cfile,'(I)')ifile
           fname = trim(adjustl(infile_no_ext))//'.'//trim(adjustl(cfile))//'.hdf5'
        else
           fname = trim(adjustl(infile_no_ext))//'.hdf5'
        endif
        
        call read_snapshot_header(isnap,ifile)
        call hdf5_open_file(snap_handle,fname,readonly=.true.)
        i = Offset(0)
        j = Offset(0)+header%NumPart_ThisFile(0)-1
        if(header%NumPart_ThisFile(0) > 0) call hdf5_read_data(snap_handle,'PartType0/ParticleIDs',ParticleIDs%Type0(i:j))
        if(header%NumPart_ThisFile(0) > 0) SnapIDs%Type0(i:j) = isnap
        i = Offset(1)
        j = Offset(1)+header%NumPart_ThisFile(1)-1
        if(header%NumPart_ThisFile(1) > 0) call hdf5_read_data(snap_handle,'PartType1/ParticleIDs',ParticleIDs%Type1(i:j))
        if(header%NumPart_ThisFile(1) > 0) SnapIDs%Type1(i:j) = isnap
        i = Offset(2)
        j = Offset(2)+header%NumPart_ThisFile(2)-1
        if(header%NumPart_ThisFile(2) > 0) call hdf5_read_data(snap_handle,'PartType2/ParticleIDs',ParticleIDs%Type2(i:j))
        if(header%NumPart_ThisFile(2) > 0) SnapIDs%Type2(i:j) = isnap
        i = Offset(3)
        j = Offset(3)+header%NumPart_ThisFile(3)-1
        if(header%NumPart_ThisFile(3) > 0) call hdf5_read_data(snap_handle,'PartType3/ParticleIDs',ParticleIDs%Type3(i:j))
        if(header%NumPart_ThisFile(3) > 0) SnapIDs%Type3(i:j) = isnap
        i = Offset(4)
        j = Offset(4)+header%NumPart_ThisFile(4)-1
        if(header%NumPart_ThisFile(4) > 0) call hdf5_read_data(snap_handle,'PartType4/ParticleIDs',ParticleIDs%Type4(i:j))
        if(header%NumPart_ThisFile(4) > 0) SnapIDs%Type4(i:j) = isnap
        i = Offset(5)
        j = Offset(5)+header%NumPart_ThisFile(5)-1
        if(header%NumPart_ThisFile(5) > 0) call hdf5_read_data(snap_handle,'PartType5/ParticleIDs',ParticleIDs%Type5(i:j))
        if(header%NumPart_ThisFile(5) > 0) SnapIDs%Type5(i:j) = isnap
        call hdf5_close_file(snap_handle)

        do itype = 0,5
           Offset(itype) = Offset(itype) + header%NumPart_ThisFile(itype)
        enddo

     enddo
  enddo

end subroutine read_hdf5_ics
       









